///\todo upgrade

//////////////////////////////////////////////////////////////////////////////////
//                                                                              //
//  This file is part of the buola project (https://code.google.com/p/buola/).  //
//                                                                              //
//  Copyright(c) 2007-2012 Xavi Gratal                                          //
//  gratal AT gmail DOT com                                                     //
//                                                                              //
//  Buola is free software: you can redistribute it and/or modify               //
//  it under the terms of the GNU General Public License as published by        //
//  the Free Software Foundation, either version 3 of the License, or           //
//  (at your option) any later version.                                         //
//                                                                              //
//  Buola is distributed in the hope that it will be useful,                    //
//  but WITHOUT ANY WARRANTY; without even the implied warranty of              //
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               //
//  GNU General Public License for more details.                                //
//                                                                              //
//  You should have received a copy of the GNU General Public License           //
//  along with buola.  If not, see <http://www.gnu.org/licenses/>.              //
//                                                                              //
//////////////////////////////////////////////////////////////////////////////////

#ifndef _BUOLA_MAT_MATRIXOPS_H_
#define _BUOLA_MAT_MATRIXOPS_H_

#include <buola/mat/cmat.h>
#include <buola/mat/initializers.h>

namespace buola { namespace mat {

template<typename tScalar>
CMat<tScalar> mpow(const CMat<tScalar> &pBase,size_t pExp)
{
    assert(pBase.Rows()==pBase.Cols());

    if(pExp==0) return mat::eye(pBase.Rows());
    if(pExp==1) return pBase;
    
    mat::CMat_d lResult;

    size_t lCurExp=1;

    if(pExp&1)
    {
        lResult=pBase;
    }
    else
    {
        lResult=pBase*pBase;
        lCurExp=2;
        
        while(!(pExp&lCurExp))
        {
            lResult*=lResult;
            lCurExp*=2;
        }
    
        if(lCurExp==pExp) return std::move(lResult);
    }
    
    mat::CMat_d lCurMat=lResult*lResult;
    lCurExp*=2;
    if(pExp&lCurExp)
        lResult*=lCurMat;
    
    while(lCurExp*2<pExp)
    {
        lCurMat*=lCurMat;
        lCurExp*=2;
        if(pExp&lCurExp)
            lResult*=lCurMat;
    }
    
    return lResult;
}
    
/*namespace mat*/ } /*namespace buola*/ }

#endif
